Numerical Exercise Nr 1

One common way to characterize the structure of a universe is to measure the statistical properties of its haloes. The simplest statistic we can perform is to count the number of haloes found for a given mass interval, normalized to the volume of the universe that contains those structures. This function is usually called the “halo mass function”. In the first part, you will compute the halo mass function for three N-body simulations (DM-only) run using a Λ-CDM cosmological scenario. You will contrast the outputs with theoretical predictions and models. For this activity you will need to install colossus using the command line: pip install colossus

1.- Donwload the data

This was uploaded to the Drive of the course and downloaded the full directory (NumEx1), the directory contains

• IllustrisTNG300-Dark at redshifts $z=$ 0, 0.5, 1, 2.   (Springel et al., 2018)
• The Small MultiDark Planck simulation at $z=$ 0.     (Klypin et al., 2016)
• The EAGLE DM-only at $z=$ 0.             (Schaye et al., 2015)

Some of this data will be show as follows

import h5py
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import numpy as np
import pandas as pd
from scipy.interpolate import interp1d 
import matplotlib.ticker as ticker

%matplotlib inline


COLOR = 'lightgray'

mpl.rcParams['text.color']      = COLOR  # labels, titles, etc colors
mpl.rcParams['axes.labelcolor'] = COLOR  # Numbers in axis colors
mpl.rcParams['xtick.color']     = COLOR  # ticks colors
mpl.rcParams['ytick.color']     = COLOR  # ticks colors
mpl.rc(      'axes',edgecolor   = COLOR) # axis Colors

mpl.rc('font', size=15) # Font Size

mpl.rcParams['figure.figsize'] = [12, 12] # Size of figures

# Function to get the values for the extend in plt.imshow
def EXTENT(x,y):
    w1=(x[1]-x[0])/2.0
    w2=(y[1]-y[0])/2.0
    return (x.min()-w1,x.max()+w1,y.min()-w2,y.max()+w2)

Read Small Multi Dark Simulation

#More info, see: https://www.cosmosim.org/metadata/mdr1/rockstar/
file_path = 'mdr1_rockstar.csv'
df = pd.read_csv(file_path)

h = 0.6777

MDR_Mvir = df[['mvir']]/h #Msol/h
MDR_X = df[['x']] #Mpc/h
MDR_Y = df[['y']] #Mpc/h
MDR_Z = df[['z']] #Mpc/h

Limit = (MDR_Mvir>1e10)
Volume= (1000/h)**3


df_filter = df[(df['z'] > 300.) & (df['z'] < 310.)]

Plot of a slice of the dark matter haloes

#plt.hist2d(df_filter['x'],df_filter['y'],bins=200,cmap='hot')
H,x,y=np.histogram2d(df_filter['x'],df_filter['y'],bins=200)
IM=plt.imshow(np.log10(H.T+1),origin='lower',extent=EXTENT(x,y),interpolation='bilinear',cmap='inferno')
CM=plt.colorbar(IM,fraction=0.046, pad=0.04,location='top')
CM.set_label(label='$Log_{10}(Count+1)$',labelpad=20)
plt.axis('scaled')
plt.xlabel('X [cMpc/h]')
plt.ylabel('Y [cMpc/h]')
print()

Read TNG300-1-Dark

redshift = [0.0, 0.1, 0.5, 1.0,  2.0] 

h     = 0.6774
Volume2= (205/h)**3 # Mpc ^3


with h5py.File('TNG300_Dark_z_'+str(redshift[0])+'.hdf5','r') as f1:
    TNG300_pos_x = f1['Pos'][:,0]
    TNG300_pos_y = f1['Pos'][:,1]
    TNG300_pos_z = f1['Pos'][:,2]
    TNG300_mvir = f1['M_Crit200'][:]/h
    Limit2 = (TNG300_mvir>1e10)

Variables=[]
for i in redshift:
    with h5py.File('TNG300_Dark_z_'+str(i)+'.hdf5','r') as f1:
        globals()['TNG300_mvir_z_'+str(i)] = f1['M_Crit200'][:]/h
        LIM = (globals()['TNG300_mvir_z_'+str(i)]>1e10)  
        globals()['TNG300_mvir_z_'+str(i)]=globals()['TNG300_mvir_z_'+str(i)][LIM]
        Variables.append('TNG300_mvir_z_'+str(i))

Plot of a slice of the dark matter haloes

fil = np.where((TNG300_mvir>1e11)&(TNG300_pos_z>100000)&(TNG300_pos_z<125000))
# plt.hist2d(TNG300_pos_x[fil]/1000.,TNG300_pos_y[fil]/1000.,bins=300,cmap='hot')
H2,x2,y2=np.histogram2d(TNG300_pos_x[fil]/1000.,TNG300_pos_y[fil]/1000.,bins=300)
IM2=plt.imshow((H2.T)**0.75,origin='lower',extent=EXTENT(x2,y2),interpolation='gaussian',cmap='inferno')
CM2=plt.colorbar(IM2,fraction=0.046, pad=0.04,location='top')
CM2.set_label(label='$(Count)^{0.75}$',labelpad=20)
plt.axis('scaled')
plt.xlabel('X [cMpc/h]')
plt.ylabel('Y [cMpc/h]')
print()

Read EAGLE 100Mpc DMO

file_path = 'EAGLE_100Mpc_DMO.txt'
df_EAGLE = pd.read_csv(file_path, delimiter=',', skiprows=20) 
EAGLE_Mvir = df_EAGLE['Group_M_Crit200']#Msol
EAGLE_X = df_EAGLE['GroupCentreOfPotential_x']#Mpc
EAGLE_Y = df_EAGLE['GroupCentreOfPotential_y']#Mpc
EAGLE_Z = df_EAGLE['GroupCentreOfPotential_z']#Mpc

Limit3 = (EAGLE_Mvir>1e10)
Volume3= 100 **3

df_EAGLE_filter = df_EAGLE[(df_EAGLE['GroupCentreOfPotential_z'] > 40.) & (df_EAGLE['GroupCentreOfPotential_z'] < 50.)]

#plt.hist2d(df_EAGLE_filter['GroupCentreOfPotential_x'],df_EAGLE_filter['GroupCentreOfPotential_y'],bins=300,cmap='hot')
H3,x3,y3=np.histogram2d(df_EAGLE_filter['GroupCentreOfPotential_x'],df_EAGLE_filter['GroupCentreOfPotential_y'],bins=300)
IM3=plt.imshow(H3.T,origin='lower',extent=EXTENT(x3,y3),interpolation='bilinear',cmap='inferno')
CM3=plt.colorbar(IM3,fraction=0.046, pad=0.04)
CM3.set_label(label='$Count$',labelpad=20)
plt.axis('scaled')
plt.xlabel('X [cMpc]')
plt.ylabel('Y [cMpc]')
print()

#==================================
# Count particules in 3 dimensions

if True: # High Resolution
    x, y, z = EAGLE_X , EAGLE_Y , EAGLE_Z
    
    N= len(x)
    #print(N)

    H, edges=np.histogramdd((x,y,z),bins=50) #bins=(50,45,55)

    x= (edges[0][1:] + edges[0][:-1])/2
    y= (edges[1][1:] + edges[1][:-1])/2
    z= (edges[2][1:] + edges[2][:-1])/2

    x,y,z = np.meshgrid(x,y,z)

    x,y,z=x.ravel(),y.ravel(),z.ravel()

    #print(len(edges[0]),len(edges[1]),len(edges[2]))
    #print(min(x),max(x),min(y),max(y),min(z),max(z))
    #print(len(x),len(y),len(z))

if 1: # Scatter density version (plotly)
    import plotly.graph_objects as go
    from plotly.subplots import make_subplots

    #------------------
    # Figure

    fig = make_subplots(rows=1, cols=1,
                    specs=[[{'type': 'surface'}]],
                    print_grid=False,
                    #subplot_titles=( f'EAGLE'),
                    )

    #------------------
    # Scatter plot High resolution

    trace1 = go.Volume(name         = 'EAGLE',
                       x            = y,
                       y            = x,
                       z            = z,
                       value        = H.ravel(),
                       #isomin      = 0.1,
                       #isomax      = 0.8,
                       opacity      = 0.1, # needs to be small to see through all surfaces
                       surface_count= 12, # needs to be a large number for good volume rendering
                       colorscale   = 'matter',
                       colorbar     = dict(thickness    = 20,
                                           outlinewidth = 0.7,
                                           title        = "Count particles per bin",
                                           titleside    = "right",
                                           tickmode     = "array",
                                           tickvals     = [], #[0.001,0.999]
                                           #labelalias   = {1: str(H.max()), 0:str(H.min()) },
                                           ticks        = "outside",
                                           x            = 1.2,)  # position of colorbar
                      )
    fig.add_trace(trace1, # trace= : figure
                  row=1,  # row=   : Position y
                  col=1)  # col=   : Position x
    
    #------------------
    # Background color (page and axes)

    # https://www.w3schools.com/colors/colors_rgb.asp

    fig['layout']['paper_bgcolor']="rgb(17,17,17)"
        
    set_bgcolor= dict(showbackground  = True,
                      backgroundcolor = "rgb(212, 219, 249)", #"rgb(101, 101, 101)"
                      gridcolor       = "rgb(10, 10, 10)",
                      zeroline        = False)

    fig.update_scenes(xaxis=set_bgcolor,
                      yaxis=set_bgcolor,
                      zaxis=set_bgcolor)
    
    #------------------
    # Font color and size

    fig.update_layout(font_color = 'white', font_size = 16)

    
    fig.update_layout(height=900, width=900, title_text="EAGLE")
    fig.show()
    print()

2.- Compute the halo mass function at $z = 0$ for the three DM-only simulations, assuming a virial mass threshold in $M_{200c} = 10^{10} M_⊙$. Compare the results obtained between the different simulations.

We have defined in the previous step the variables that contain the mass of the halos for each simulation and we also make the conversion to obtain this value without solar masses when applicable.

# Data Counted # This step takes a while

#print(sum(Limit.to_numpy())[0],(MDR_Mvir.to_numpy())[:,0].shape)
h1,x1=np.histogram(np.log10(MDR_Mvir[Limit]),bins=100)           # We count the data in 100 bins log spaced
w1= (x1[1]-x1[0])                                                # Width of the bins (mass enclosed)
X1= (x1[1:]+x1[:-1])/2.                                          # Center of the bins

#print(sum(Limit2),TNG300_mvir.shape)
h2,x2=np.histogram(np.log10(TNG300_mvir[Limit2]),bins=100)
w2= (x2[1]-x2[0])
X2= (x2[1:]+x2[:-1])/2.

#print(sum(Limit3),EAGLE_Mvir.shape)
h3,x3=np.histogram(np.log10(EAGLE_Mvir[Limit3]),bins=100)
w3= (x3[1]-x3[0])
X3= (x3[1:]+x3[:-1])/2.

We show just the count with the impossed limit

plt.style.use('dark_background')
plt.style.use('seaborn-bright')

plt.stairs(h1,x1,fill=False,label='Small MultiDark Planck simulation')

plt.stairs(h2,x2,fill=False,label='IllustrisTNG300-Dark')

plt.stairs(h3,x3,fill=False,label='EAGLE DM-ONLY')

plt.yscale('log')
plt.ylabel('$N_{halo}$')
plt.xlabel('$\log_{10}(M [M_\odot])$')
plt.legend()
plt.grid(ls=":",color='lightgray',lw=0.5,which='both')

We now show the halo mass function

(The variables of volume were defined when we read the data)

plt.plot(X1,h1/w1/Volume,marker='o',label='Small MultiDark Planck simulation')

plt.plot(X2,h2/w2/Volume2,marker='d',label='IllustrisTNG300-Dark')

plt.plot(X3,h3/w3/Volume3,marker='s',label='EAGLE DM-ONLY')

plt.yscale('log')
#plt.xscale('log')
plt.ylabel(r'$N_{halo}\cdot\ dv^{-1}\cdot\ d(Log_{10}(M_{200}))^{-1}$',fontsize=20)
plt.xlabel('$\log_{10}(M [M_\odot])$')

plt.legend(loc='upper left', bbox_to_anchor=(1.05, 1.))
plt.grid(ls=":",color='lightgray',lw=0.5,which='both')

3.- Usecollossus (https://bdiemer.bitbucket.io/colossus/lss_mass_function.html) to compare the results with the theoretical model from the assumption of spherical and ellipsoidal collapse (Press & Schechter 1974, Sheth & Tormen 1999).

from colossus.cosmology import cosmology
from colossus.lss import mass_function

we will create the distributions and make an interpolation of them in order to make much easier the comparison

cosmology.setCosmology('WMAP9')

new_x   = np.linspace(9.99,16.01,100)

mfunc_P = mass_function.massFunction(10**new_x, 0.0, mdef = 'fof', model = 'press74',q_out='dndlnM')
mfunc_S = mass_function.massFunction(10**new_x, 0.0, mdef = 'fof', model = 'sheth99',q_out='dndlnM')
# EXTRA   = mass_function.massFunction(10**new_x, 0.0, mdef = 'vir', model = 'tinker08',q_out='dndlnM')

MH1= interp1d(new_x,mfunc_P)
MH2= interp1d(new_x,mfunc_S)

Y1=h1/w1/Volume

Y1[Y1<=0]=np.nan # easier to handle

Y2=h2/w2/Volume2

Y2[Y2<=0]=np.nan # easier to handle

Y3=h3/w3/Volume3

Y3[Y3<=0]=np.nan # easier to handle

Just Ploting them directly

plt.plot(X1,Y1,lw=3,label='Small MultiDark Planck simulation',color='b')#,marker='o')

plt.plot(X2,Y2,lw=3,label='IllustrisTNG300-Dark',color='springgreen')#,marker='d')

plt.plot(X3,Y3,lw=3,label='EAGLE DM-ONLY',color='r')#,marker='s')



plt.plot(new_x,mfunc_P,ls='--',lw=3,label='Press and Schechter',color='fuchsia')#,marker='.')

plt.plot(new_x,mfunc_S,ls='--',lw=3,label='Sheth and Tormen',color='aqua')#,marker='*')

# plt.plot(new_x,EXTRA,label='EXTRA',marker='p')



plt.ylabel(r'$N_{halo}\cdot\ dv^{-1}\cdot\ d(Log_{10}(M_{200}))^{-1}$',fontsize=20)
plt.xlabel('$\log_{10}(M [M_\odot])$')


plt.legend(loc='upper left', bbox_to_anchor=(1.02, 1.))
plt.grid(ls=":",color='lightgray',lw=0.5,which='both')
plt.ylim(2e-8,2)
plt.yscale('log')

But this is hard to see since is very crowded so let's see it in another way

fig= plt.figure(figsize=(18,9))
gs = gridspec.GridSpec( 3,3,                            # Number of axis y,x
                        height_ratios = [0.5,0.05,0.5],  # relatives ratios of heigh
                        width_ratios  = [0.33,0.33,0.33],               # relatives ratio of with
                        left  = 0.1,      # Space to the edge of left   from the nearest axis
                        right = 0.95,     # Space to the edge of right  from the nearest axis
                        bottom= 0.1,     # Space to the edge of bottom from the nearest axis
                        top   = 0.97,     # Space to the edge of top    from the nearest axis
                        wspace= 0.15,      # Space horizontal between each of the axis
                        hspace= 0.0)      # Space vertical   between each of the axis

ax1  = fig.add_subplot(gs[0,0])
ax2  = fig.add_subplot(gs[2,0],sharex=ax1,sharey=ax1)
ax3  = fig.add_subplot(gs[0,1],sharex=ax1,sharey=ax1)
ax4  = fig.add_subplot(gs[2,1],sharex=ax1,sharey=ax1)
ax5  = fig.add_subplot(gs[0,2],sharex=ax1,sharey=ax1)
ax6  = fig.add_subplot(gs[2,2],sharex=ax1,sharey=ax1)

# (here we can define if we want to share axes and with which one)
#--
# Ploting

Center=1.0

Op= lambda a,b: a/b

ax1.set_yscale('log')
ax1.yaxis.set_major_formatter(ticker.FuncFormatter(lambda y,pos: ('{{:.{:1d}f}}'.format(int(np.maximum(-np.log10(y),0)))).format(y)))

fig.supxlabel('$\log_{10}(M [M_\odot])$',fontsize=18)
fig.supylabel('Ratio',fontsize=18)
#fig.supylabel('Residual',fontsize=18)


# Press and Schechter
ax1.plot(X1,Op(Y1,MH1(X1)),lw=3,label='Small MultiDark Planck simulation',color='b')
ax1.axhline(Center,ls='--',lw=3,label='Press and Schechter',color='fuchsia',alpha=0.5)
ax1.plot([],[],ls='--',lw=3,label='Sheth and Tormen',color='aqua',alpha=0.5)
ax1.legend(loc='upper center', bbox_to_anchor=(.5, 1.25))


ax3.plot(X2,Op(Y2,MH1(X2)),lw=3,label='IllustrisTNG300-Dark',color='springgreen')
ax3.axhline(Center,ls='--',lw=3,label='Press and Schechter',color='fuchsia',alpha=0.5)
ax3.plot([],[],ls='--',lw=3,label='Sheth and Tormen',color='aqua',alpha=0.5)
ax3.legend(loc='upper center', bbox_to_anchor=(.5, 1.25))

ax5.plot(X3,Op(Y3,MH1(X3)),lw=3,label='EAGLE DM-ONLY',color='r')
ax5.axhline(Center,ls='--',lw=3,label='Press and Schechter',color='fuchsia',alpha=0.5)
ax5.plot([],[],ls='--',lw=3,label='Sheth and Tormen',color='aqua',alpha=0.5)
ax5.legend(loc='upper center', bbox_to_anchor=(.5, 1.25))


# Sheth and Tormen
ax2.plot(X1,Op(Y1,MH2(X1)),lw=3,label='Small MultiDark Planck simulation',color='b')
ax2.axhline(Center,ls='--',lw=3,label='Sheth and Tormen',color='aqua',alpha=0.5)
# ax2.legend(loc='upper right')

ax4.plot(X2,Op(Y2,MH2(X2)),lw=3,label='IllustrisTNG300-Dark',color='springgreen')
ax4.axhline(Center,ls='--',lw=3,label='Sheth and Tormen',color='aqua',alpha=0.5)
# ax4.legend(loc='upper right')

ax6.plot(X3,Op(Y3,MH2(X3)),lw=3,label='EAGLE DM-ONLY',color='r')
ax6.axhline(Center,ls='--',lw=3,label='Sheth and Tormen',color='aqua',alpha=0.5)
# ax6.legend(loc='upper right')



# Erasing extra y labels
plt.setp(ax1.get_xticklabels(), visible=False)
plt.setp(ax3.get_xticklabels(), visible=False)
plt.setp(ax5.get_xticklabels(), visible=False)
print()

4.- Interpret the results and differences between the simulations and theoretical models. Con- sider the differences in mass resolution, and size of the simulated box.

From this last plot we can clearly see that the ellipsoidal collapse fits better over all the simulation, now if we look more in detail:

• The SMP particulary gets better fits at larges masses from the ellipsoidal collapse while we can see that for lower masses non of them is able to reproduce the behaviour and this makes sense since this it's likely a problem of the resolution of the box since is the largest of our data set

• For illustris TNG300 we have that the ellipsoidal collapse represents very well the trend of the data even if at the largest masses tends to deviate this could be easily associated with the sise of the box

• While for EAGLE the data is very well represent by an ellipsoidal collapse

5.- For IllustrisTNG300-DMO, compute the halo mass function in the four redshifts provided. Fit the halo mass function using the following expression:

$$\frac{dN}{dM}\sim\left(\frac{M_{halo}}{m_0}\right)^n$$

measure the slope n of the halo mass function. How does the slope vary as a function of redshift? What can you conclude about the way that galaxies build up over time ?

# We read the variables (they were defined when we read the TNG-300-dark-1)
Sec_Var={}

for i in Variables:
    D=globals()[i]
    #print(D.shape)
    HH,XX=np.histogram(np.log10(D),bins=100)
    mask= (HH>0)
    
    WW   = (XX[1]-XX[0])
    XX   = ((XX[1:]+XX[:-1])/2.)[mask]
    
    Sec_Var[i]=np.array([XX,HH[mask]/WW/Volume2])
    

fig= plt.figure(figsize=(18,9))
gs = gridspec.GridSpec( 1,3,                            # Number of axis y,x
                        height_ratios = [1],  # relatives ratios of heigh
                        width_ratios  = [0.5,0.1,0.5],               # relatives ratio of with
                        left  = 0.1,      # Space to the edge of left   from the nearest axis
                        right = 0.95,     # Space to the edge of right  from the nearest axis
                        bottom= 0.1,     # Space to the edge of bottom from the nearest axis
                        top   = 0.97,     # Space to the edge of top    from the nearest axis
                        wspace= 0.15,      # Space horizontal between each of the axis
                        hspace= 0.0)      # Space vertical   between each of the axis

ax1  = fig.add_subplot(gs[0,0])
ax2  = fig.add_subplot(gs[0,2],sharex=ax1,sharey=ax1)

print('File\t\t\t[Slope (N)    Coeff (M0)]')
print('-------------------------------------------------')
for i in Sec_Var.keys():
    D = Sec_Var[i]
    ax1.plot(D[0],D[1])
    
    pf=np.polyfit(D[0],np.log10(D[1]),1)
    
    print(f'{i}\t{pf}')
    

    ax2.plot(D[0],10**((D[0]*pf[0])+pf[1]),label=i)

ax1.set_yscale('log')
ax2.legend()

ax1.grid(zorder=0,ls=':',color='gray')
ax2.grid(zorder=0,ls=':',color='gray')

fig.supylabel(r'$N_{halo}\cdot\ dv^{-1}\cdot\ d(Log_{10}(M_{200}))^{-1}$',fontsize=20)
fig.supxlabel('$\log_{10}(M [M_\odot])$',fontsize=18)

print()

File			[Slope (N)    Coeff (M0)]
-------------------------------------------------
TNG300_mvir_z_0.0	[-1.01584723  9.7721955 ]
TNG300_mvir_z_0.1	[-1.01394106  9.75779471]
TNG300_mvir_z_0.5	[-1.06472794 10.3614411 ]
TNG300_mvir_z_1.0	[-1.12013746 10.98513452]
TNG300_mvir_z_2.0	[-1.24970781 12.3532928 ]

The slope presents slighty changues over redshift , in specific this decrease with redshift and this makes sense if we assume the merger trees evolution of the universe (we have less massive halos at early stages of the universe and more of small ones)

6.- Using the snapshots from IllustrisTNG300-DMO, compute the number density of dark matter haloes for three halo mass bins, ($10^{10}$ − $10^{11}$, $10^{11}$ − $10^{12}$, $10^{12}$ − $10^{13}$) $M_⊙$ as a function of redshift. Interpret the results in terms of the building up of structures within the cosmological paradigm adopted.

masses = [1e10, 1e11, 1e12, 1e13]

Data = []

for d in Variables: # Getting the info of the specific redshift 
    D    = globals()[d] 
    aux  = []
    for i in range(len(masses)-1):
        l,u  = masses[i],masses[i+1] # Range of masses desired 
        mask = (l<D)&(D<u)
        aux.append(len(D[mask])/Volume2/(np.log10(u)-np.log10(l))) # Number of halos/ Volume / bin mass width
        
    Data.append(aux)

Data=np.array(Data)

fig= plt.figure(figsize=(18,6))
gs = gridspec.GridSpec( 1,5,                             # Number of axis y,x
                        height_ratios = [1],             # relatives ratios of heigh
                        width_ratios  = [1,0.1,1,0.1,1], # relatives ratio of with
                        left  = 0.1,      # Space to the edge of left   from the nearest axis
                        right = 0.95,     # Space to the edge of right  from the nearest axis
                        bottom= 0.15,     # Space to the edge of bottom from the nearest axis
                        top   = 0.97,     # Space to the edge of top    from the nearest axis
                        wspace= 0.1,      # Space horizontal between each of the axis
                        hspace= 0.0)      # Space vertical   between each of the axis

ax1  = fig.add_subplot(gs[0,0])
ax2  = fig.add_subplot(gs[0,2])
ax3  = fig.add_subplot(gs[0,4])

AX=[ax1,ax2,ax3]


# Function to put scientific notation
def nsf(x):
    return f'10^{{{np.log10(x):.0f}}}'

# Ploting the masses
for i in range(0,len(masses)-1):
    l,u  = masses[i],masses[i+1] # Range of masses desired 
    DDD  = Data[:,i]
    AX[i].plot(redshift,DDD,marker='o',label=f'${nsf(l)}-{nsf(u)}$')

    AX[i].legend(loc='upper left', bbox_to_anchor=(0.52, 1.12))
    AX[i].ticklabel_format(axis='y', scilimits=[-1, 1])
    AX[i].grid(zorder=0,ls=':',color='gray')
    #AX[i].set_yscale('log')

fig.supylabel(r'$N_{halo}\cdot\ dv^{-1}\cdot\ d(Log_{10}(M_{200}))^{-1}$',fontsize=20)
fig.supxlabel('Redshift',fontsize=20)
print()

This behaviour can be easily understood in the context of the hierarchical universe, this means that the amount of small halos (first panel) will naturally decreases as function of redshift due they merge along the time , in the other extreme the most massive ones (last panel) will be more abundant in later stages of the universe since more halos have had enough time to being formed, finally the middle panel reach a peak at certain redshift where they have just enough time to being formed but not much to be merged.

print('Datapoint in the last graph')
Data

Datapoint in the last graph

array([[0.10847032, 0.01450199, 0.001861  ],
       [0.11198504, 0.0149035 , 0.00189412],
       [0.12394938, 0.0161267 , 0.00197   ],
       [0.13284127, 0.0167439 , 0.00188525],
       [0.13456275, 0.01507434, 0.00128826]])